################################################################################################
###Replication Code for "Because He's Gay": How Race, Gender, and Sexuality Shape Perceptions of Judicial Fairness
###Journal of Politics
###Authors: Ana Bracic, Mackenzie Israel-Trummel, Tyler Johnson, and Kathleen Tipler
################################


################################
##Instructions:
##This code file is to analyze Democratic respondents.
##First run the code file 1_Conjoint Set Up to generate the stacked data file for analysis.
################################

##Restrict to Democrats. Load stacked data. 
cd <- read.csv(file = "stackedconjointdata.csv")
length(cd$id) #40420 observations total
table(cd$Rpid)

dems <- subset(cd, cd$Rpid < 4)
table(dems$Rpid)
length(dems$id) #16,500 observations of Dems
names(dems)

dems$race <- as.factor(dems$race)

dems$race <- relevel(dems$race, ref="White")


dems$sexuality <- as.factor(dems$sexuality)
dems$sexuality <- relevel(dems$sexuality, ref="Straight")


dems$lawschool <- as.factor(dems$lawschool)
dems$lawschool <- factor(dems$lawschool, levels = c("Top 10, Public", "Top 10, Private", "#11-25, Private", "#11-25, Public", "#26-50, Private" ,   "#26-50, Public",  "#51-100, Private",   "#51-100, Public", "#101-200, Private",  "#101-200, Public"))

dems$occupation <- as.factor(dems$occupation)
dems$occupation <-  relevel(dems$occupation, ref="Law professor")

dems$Rrace <- as.factor(dems$Rrace)
dems$Rrace <- relevel(dems$Rrace, ref="White")



## Load necessary packages for importing the function for robust clustered errors
library(RCurl)

## Import the function from repository
url_robust <- "https://raw.githubusercontent.com/IsidoreBeautrelet/economictheoryblog/master/robust_summary.R"
eval(parse(text = getURL(url_robust, ssl.verifypeer = FALSE)),
     envir=.GlobalEnv)


## New function: summary(modelname, cluster = c("id)) provides you with robust clustered error adjustment


##############################
### Regressions and Effects Figures
###############################

########################################################################
## CREATE FIGURE 2A
########################################################################

## Conservatism of judge
modC.D <- lm(conservative_1 ~ race + gender + sexuality + party + agecat + lawschool + occupation
            + Rwoman + RLGB + Rincome + Rrelig + Rlibcon + Rrace + Redu + Rage + Rknow    
            +factor(vignette_id),
            data = dems)

summary(modC.D)

## use new summary function
summary(modC.D)
modC2.D <- summary(modC.D, cluster = c("id"))



coefsC.D <- c(modC2.D$coef[,1]) #Pull coefficients from model
confintC.D <- c(modC2.D$coef[,2]*1.96) #Create 95 percent confidence interval from standard error 

## Join into data frame
tmpC.D <- data.frame(cbind(coefsC.D,confintC.D))



##Create slimmed data frame to just attributes for plot

tmpC.D2 <- tmpC.D[grep("^race|^genderWoman$|^sexualityGay$|^partyRepublican",row.names(tmpC.D)),]
tmpC.D2$variable <- c("Asian", "Black", "Hispanic", "Native", "Woman", "Gay", "Republican")
tmpC.D2$variable <- as.factor(tmpC.D2$variable)
names(tmpC.D2) <- c("coef", "confint", "variable")


## Load list of variables for plot axis
variablesC.D <- read.csv(file = "conjointvariablesSHORT.csv")


## Input data for plotting
for (i in 1:dim(variablesC.D)[1]){
  if (sum(tmpC.D2$variable == variablesC.D$variable[i], na.rm = T) > 0){
    variablesC.D$coef[i] <- tmpC.D2$coef[tmpC.D2$variable == variablesC.D$variable[i]]
    variablesC.D$confint[i] <- tmpC.D2$confint[tmpC.D2$variable == variablesC.D$variable[i]]
  }
}


## Reorder for plotting
variables2C.D <- variablesC.D

order <- 1:dim(variables2C.D)[1]
variables2C.D <- transform(variables2C.D, variable2 = reorder(variable2, -order) ) 

variables2C.D$labels <- round(variables2C.D$coef, 3)
variables2C.D$labels[variables2C.D$labels == 0] <- NA

##PLOT

plotC.D <- ggplot(data = variables2C.D, aes(x = coef, y = variable2)) +
  geom_point(aes(colour = group, shape = group), size = 4.5)+
  geom_text(aes(label = labels), nudge_y = 0.4, nudge_x = -0.002, size=6)+
  geom_errorbarh( aes(y = variable2, xmin = coef - confint, xmax = coef + confint, height=.28,colour=group),size=.5)+
  scale_color_manual(values=c("#376795", "#376795", "#376795", "#376795")) + 
  scale_shape_manual(values=c(16, 16, 16,  16)) + 
  theme_bw() + 
  xlim(-0.1, 0.2) + 
  xlab("Change in Perceived Conservatism")+
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_text(size = 20)) +
  geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") +
  theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black", size = 20))  +
  theme(axis.text.x = element_text(size = 20))  +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(0.1, .5, 0.1, .5), "cm"))


pdf(file = "Fig2A.pdf", height=9, width=8)
plotC.D
dev.off()




########################################################################
## CREATE FIGURE 3A
########################################################################

## Selected to hear case fairly


modS.D <- lm(selected ~ race + gender + sexuality + party + agecat + lawschool + occupation
            + Rwoman + RLGB + Rincome + Rrelig + Rlibcon + Rrace + Redu + Rage + Rknow    
            + factor(vignette_id),
            data = dems)

# use new summary function
summary(modS.D)
modS2.D <- summary(modS.D, cluster = c("id"))


coefsS.D <- c(modS2.D$coef[,1]) #Pull coefficients
confintS.D <- c(modS2.D$coef[,2]*1.96) #Create 95 percent confidence interval from standard error 

## Join into data frame
tmpS.D <- data.frame(cbind(coefsS.D, confintS.D))

## Create slimmed data frame to just attributes for plot
tmpS.D2 <- tmpS.D[grep("^race|^genderWoman$|^sexualityGay$|^partyRepublican",row.names(tmpS.D)),]
tmpS.D2$variable <- c("Asian", "Black", "Hispanic", "Native", "Woman", "Gay", "Republican")
tmpS.D2$variable <- as.factor(tmpS.D2$variable)
names(tmpS.D2) <- c("coef", "confint", "variable")



#input list of variables for plot axis
variablesS.D <- read.csv(file = "conjointvariablesSHORT.csv")

#input data for plotting
for (i in 1:dim(variablesS.D)[1]){
  if (sum(tmpS.D2$variable == variablesS.D$variable[i], na.rm=T)>0){
    variablesS.D$coef[i] <- tmpS.D2$coef[tmpS.D2$variable == variablesS.D$variable[i]]
    variablesS.D$confint[i] <- tmpS.D2$confint[tmpS.D2$variable == variablesS.D$variable[i]]
  }
}

#reorder for plotting
variables2S.D <- variablesS.D


order <- 1:dim(variables2S.D)[1]
variables2S.D <- transform(variables2S.D, variable2 = reorder(variable2, -order) ) 

table(variables2S.D$variable2 , variables2S.D$group)
variables2S.D$labels <- formatC(variables2S.D$coef, format = "f", digits = 3)
variables2S.D$labels[variables2S.D$labels == "0.000"] <- NA


#PLOT

plotS.D <- ggplot(data = variables2S.D, aes(x = coef, y = variable2)) +
  geom_point(aes(shape = group, color = group), size = 4.5)+
  geom_text(aes(label = labels), nudge_y = 0.35, size=6)+
  geom_errorbarh( aes(y = variable2, xmin = coef - confint, xmax = coef + confint, height=.28,colour = group),size = .5)+
  scale_color_manual(values = c("#376795", "#376795", "#376795", "#376795")) + 
  scale_shape_manual(values = c(16, 16, 16,  16)) + 
  theme_bw() + 
  xlab("Change in Judge Preference")+
  xlim(-.17, .15)+
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_text(size = 20)) +
  geom_vline(xintercept = 0,size=.5,colour = "black",linetype = "dotted") +
  theme(axis.text.y = element_text(angle = 0, hjust = 0, color = "black", size = 20))  +
  theme(axis.text.x = element_text(size = 20))  +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(0.1, .5, 0.1, .5), "cm"))


pdf(file = "Fig3A.pdf", height=9, width=8)
plotS.D
dev.off()


########################################################################
## CREATE FIGURE 4A
########################################################################

## Politically motivated model

modP.D <- lm(political_1 ~ race + gender + sexuality + party + agecat + lawschool + occupation
            + Rwoman + RLGB + Rincome + Rrelig + Rlibcon + Rrace + Redu + Rage + Rknow    
            +factor(vignette_id),
            data = dems)

summary(modP.D)

# use new summary function
summary(modP.D)
modP2.D <- summary(modP.D, cluster = c("id"))


coefsP.D <- c(modP2.D$coef[,1]) #Pull coefficients
confintP.D <- c(modP2.D$coef[,2]*1.96) #Create 95 percent confidence interval from standard error 

## Join into data frame
tmpP.D <- data.frame(cbind(coefsP.D, confintP.D))

## select particular coefficients
tmpP.D2 <- tmpP.D[grep("^race|^genderWoman$|^sexualityGay$|^partyRepublican", row.names(tmpP.D)),]
tmpP.D2$variable <- c("Asian", "Black", "Hispanic", "Native", "Woman", "Gay", "Republican")
tmpP.D2$variable <- as.factor(tmpP.D2$variable)
names(tmpP.D2) <- c("coef", "confint", "variable")



#input list of variables for plot axis
variablesP.D <- read.csv(file = "conjointvariablesSHORT.csv")


#input data for plotting
for (i in 1:dim(variablesP.D)[1]){
  if (sum(tmpP.D2$variable== variablesP.D$variable[i], na.rm=T)>0){
    variablesP.D$coef[i]<- tmpP.D2$coef[tmpP.D2$variable== variablesP.D$variable[i]]
    variablesP.D$confint[i]<- tmpP.D2$confint[tmpP.D2$variable==variablesP.D$variable[i]]
  }
}

#reorder for plotting
variables2P.D <- variablesP.D

order <- 1:dim(variables2P.D)[1]
variables2P.D <- transform(variables2P.D, variable2 = reorder(variable2, -order) ) 

variables2P.D$labels <- formatC(variables2P.D$coef, format = "f", digits = 3)
variables2P.D$labels[variables2P.D$labels == "0.000"] <- NA

#PLOT

plotP.D <- ggplot(data = variables2P.D, aes(x = coef, y = variable2)) +
  geom_point(aes(colour = group, shape = group), size = 4.5)+
  geom_text(aes(label = labels), nudge_y = 0.35, size=6)+
  geom_errorbarh( aes(y = variable2, xmin = coef - confint, xmax = coef + confint, height=.28, colour=group), size=.5)+
  scale_color_manual(values=c("#376795", "#376795", "#376795", "#376795")) + 
  scale_shape_manual(values=c(16, 16, 16,  16)) + 
  theme_bw() + 
  xlab("Change in Perceived Political Motivation")+
  xlim(-0.06, 0.06)+
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_text(size = 20)) +
  geom_vline(xintercept = 0, size=.5, colour="black", linetype="dotted") +
  theme(axis.text.y = element_text(angle = 0, hjust = 0, color="black", size = 20))  +
  theme(axis.text.x = element_text(size = 20))  +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(0.1, .5, 0.1, .5), "cm"))


pdf(file = "Fig4A.pdf", height = 9, width = 8)
plotP.D
dev.off()



########################################################################
## CREATE FIGURE 5A
########################################################################



## Trust in supreme court model
modT.D <- lm(trustcourt_1 ~ race + gender + sexuality + party + agecat + lawschool + occupation
            + Rwoman + RLGB + Rincome + Rrelig + Rlibcon + Rrace + Redu + Rage + Rknow    
            +factor(vignette_id),
            data = dems)

summary(modT.D)

# use new summary function
modT2.D <- summary(modT.D, cluster = c("id"))

coefsT.D <- c(modT2.D$coef[,1])
confintT.D <- c(modT2.D$coef[,2]*1.96)

tmpT.D <- data.frame(cbind(coefsT.D, confintT.D))

##Create slimmed data frame to just attributes for plot

tmpT.D2 <- tmpT.D[grep("^race|^genderWoman$|^sexualityGay$|^partyRepublican", row.names(tmpT.D)),]
tmpT.D2$variable <- c("Asian", "Black", "Hispanic", "Native", "Woman", "Gay", "Republican")
tmpT.D2$variable <- as.factor(tmpT.D2$variable)
names(tmpT.D2) <- c("coef", "confint", "variable")


#input list of variables for plot axis
variablesT.D <- read.csv(file = "conjointvariablesSHORT.csv")


#input data for plotting
for (i in 1:dim(variablesT.D)[1]){
  if (sum(tmpT.D2$variable == variablesT.D$variable[i], na.rm=T)>0){
    variablesT.D$coef[i] <- tmpT.D2$coef[tmpT.D2$variable == variablesT.D$variable[i]]
    variablesT.D$confint[i] <- tmpT.D2$confint[tmpT.D2$variable == variablesT.D$variable[i]]
  }
}


#reorder for plotting

variables2T.D <- variablesT.D

order <- 1:dim(variables2T.D)[1]
variables2T.D <- transform(variables2T.D, variable2 = reorder(variable2, -order) ) 

variables2T.D$labels <- round(variables2T.D$coef, 3)
variables2T.D$labels[variables2T.D$labels == 0] <- NA

#PLOT

plotT.D <- ggplot(data = variables2T.D, aes(x = coef, y = variable2)) +
  geom_point(aes(colour = group, shape = group), size = 4.5)+
  geom_text(aes(label = labels), nudge_y = 0.35, size = 6) +
  geom_errorbarh(aes(y = variable2, xmin = coef - confint, xmax = coef + confint, height = .28,colour = group),size = .5)+
  scale_color_manual(values = c("#376795", "#376795", "#376795", "#376795")) + 
  scale_shape_manual(values = c(16, 16, 16,  16)) + 
  theme_bw() + 
  xlab("Change in Trust") +
  xlim(-.1, .1) + 
  theme(axis.title.y = element_blank()) +
  theme(axis.title.x = element_text(size = 20)) +
  geom_vline(xintercept = 0, size = .5, colour = "black", linetype = "dotted") +
  theme(axis.text.y = element_text(angle = 0, hjust = 0, color = "black", size = 20))  +
  theme(axis.text.x = element_text(size = 20))  +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(0.1, .5, 0.1, .5), "cm"))


pdf(file = "Fig5A.pdf", height = 9, width = 8)
plotT.D
dev.off()



